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ABSTRACT 

The  optimization  of  spacecraft  trajectories  in  vacuum  has 
received  extensive  consideration  since  the  inception  of  space 
flight,  yet,  the  effects  of  atmosphere  have  been  largely 
neglected.  The  advent  of  low  Earth-orbiting,  large  satellites 
and  platforms  necessitates  that  atmosphere  be  included  in  the 
optimization  process.  A  practical  means  of  studying  this  topic 
is  as  a  problem  in  minimum- fuel  orbital  maintenance.  Optimal 
control  theory  advances  the  notion  that  orbital  maintenance  is 
optimized  through  periodic  thrusting  as  opposed  to  forcing 
Keplarian  motion  by  nullifying  the  effects  of  drag  with 
thrust.  Further,  this  must  be  optimized  by  primer  vectoring. 
This  thesis  examined  the  efficiency  of  a  simple  method  of 
orbital  maintenance  using  fixed-angle  transverse  thrusting. 
Results  show  that  for  the  purpose  of  fuel-minimization,  the 
width  of  the  radial  band  in  which  the  satellite  is  to  be 
maintained,  is  dependent  upon  thruster  size.  In  nearly  all 
cases,  a  thrust-angle  of  70  degrees  maximized  the  fuel  saved. 
This  thesis  shows  that  fixed-angle  transverse  thrusting  does 
not  improve  on  forced  Keplarian  motion  and  hence  thrust 
vectoring  must  be  optimized. 
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I .  INTRODUCTION 

In  1963,  D.F.  Lawden,  in  his  monograph  [Ref.  1], 
laid  the  foundation  for  what  has  become  a  highly  sophisticated 
subdiscipline  of  astrodynamics,  optimization  of  space 
trajectories.  In  1979,  Marec  [Ref.  2]  provided  a  more 
comprehensive  treatment  of  the  subject  in  his  text  on  optimal 
space  trajectories.  Examination  of  the  optimization  of 
spacecraft  trajectories  has  been  treated  by  many  authors  in 
manners  similar  to  these  two  great  works,  yet,  until  recently, 
the  study  of  atmospheric  effects  upon  these  trajectories  was 
largely  neglected.  Research  and  development  of  hypervelocity 
vehicles  have  kindled  interest  in  this  area,  through  which 
study,  other  areas  of  interest  have  emerged. 

One  such  area  is  the  effect  of  aerodynamic  force  on  non- 
lifting,  or,  blunt  bodies.  First  addressed  by  Ross  and  Melton 
[Ref.  3],  this  subject  is  of  particular  interest  for 
two  reasons.  First,  atmospheric  effects  on  low-Earth  orbiting 
(LEO)  satellites  are  of  obvious  interest  and  second,  as  stated 
in  their  paper  [Ref.  3:p.  2],  better  understanding  of  this 
phenomenon  could  provide  deeper  insight  into  the  more 
complicated  topic  of  lifting  bodies  in  the  upper  atmosphere. 

Much  work  has  been  accomplished  in  this  area  pertaining  to 
accurate  prediction  of  satellite  orbits  [Ref.  4,  5] . 
The  focus  on  atmospheric  effects  as  they  pertain  to  the 


specific  problem  of  minimum- fuel  orbital  transfers,  however, 
is  unique  to  [Ref.  3]  and  the  follow-on  work  described  here. 
The  particular  problem  of  orbital  maintenance  can  provide 
insight  into  the  more  general,  and  complex  problem  of  orbital 
transfer. 

Historically,  orbital  transfers  (coplanar)  have  been 
accomplished  as  either,  two,  or,  three-impulse  maneuvers 
[Ref.  6:pp  78-88].  The  problem  of  orbital  transfer  is 
approached  as  a  minimization  of  energy  required  to  move  a 
satellite  from  orbit  A  to  orbit  B,  or,  equivalently,  a 
minimization  of  the  characteristic  velocity.  A  satellite  that 
has  descended  from  an  initial  orbit  due  to  a  disturbing  force 
such  as  drag,  and  which  must  be  returned  to  its  initial  orbit 
can  be  approached  as  a  problem  requiring  an  orbital  transfer. 
Orbital  transfers  such  as  this  are  optimally  performed  by  two- 
impulse  transfers,  such  as  a  Hohman  transfer.  Large  orbital 
transfers  (ra>11.8rp)  are  optimized  with  three-impulse 
maneuvers  [Ref.  6:p  87].  As  low-Earth  orbits  become  more 
frequently  utilized,  deeper  understanding  of  the  effects  of 
drag  must  be  achieved  in  order  that  orbits,  propulsion  systems 
and  costs  are  optimized.  This  is  particularly  applicable  for 
large  satellites,  such  as  the  proposed  space  station,  that 
must  remain  in  low  orbits  for  extended  periods  of  time. 

Additionally,  more  complex  areas  of  study,  such  as  that  of 
lifting  bodies  in  the  upper  atmosphere,  could  benefit  from  the 


insights  gained  through  a  deeper  understanding  of  atmospheric 
effects  on  non-lifting  bodies. 

Research  on  the  atmospheric  effects  on  low-Earth  orbiting 
spacecraft  is  sure  to  receive  much  attention  in  the  future. 
The  benefits  to  existing  and  future  systems,  while  extensive, 
remain  relatively  unexplored  and  demand  the  attention  of  the 
aerospace  industry. 


II.    GENERAL  FORMULATION  OF  THE  PROBLEM 

Optimization  of  orbital  transfers  is  a  subject  that  has 
achieved  a  high  degree  of  sophistication  and  many  elegant 
solutions  exist  [Ref.  1,  2].  However,  the  specific  treatment 
of  non-lifting  bodies  is  in  the  initial  stages  of  development. 
The  increasing  number  of  low  Earth  orbiting  satellites 
requires  that  a  study  of  atmospheric  effects  on  orbital 
trajectories  be  conducted.  In  this  thesis,  the  problem  of 
minimum  fuel  orbital  maintenance  is  considered.  Two  methods 
are  examined  by  which  orbital  maintenance  may  be  performed. 

One  method  is  to  counter  drag  with  thrust.  In  this  forced 
Keplarian  motion,  the  reaction  control  system  would  thrust 
continuously  for  the  duration  of  the  satellite's  lifetime  with 
magnitude  and  direction  equal  and  opposite  to  drag.  The  second 
scheme  considered  here,  utilizes  periodic  transverse 
thrusting,  or,  non-Keplarian  motion  to  correct  for 
perturbations  due  to  drag.  The  question  to  be  answered  is 
whether  or  not  optimal  non-Keplarian  trajectories  are  superior 
to  f orced-Keplarian  trajectories. 

Let  the  problem  be  defined  as  maintaining  an  orbital 
deviation  within  a  specified  radial  band  rmir<r<rmax.  Is  forced- 
Keplarian  orbital  maintenance,  i.e.,  exactly  countering 
aerodynamic  forces  with  thrust,  superior  to  non-Keplarian 
orbital   maintenance,  i.e.,  allowing  the  orbit  to  decay  to  a 


certain  point,  then  reboosting  to  a  point  above  the  desired 
altitude.  While  maintaining  the  spacecraft  within  the 
specified  radial  band?  The  study  performed  by  Ross  and  Melton 
[Ref.  3:p.  4]  suggests  that  f orced-Keplarian  motion  is  not 
optimal  and  that  thrust  vectoring  must  be  considered  if  an 
optimal  solution  is  to  be  obtained. 

This  thesis  addresses  an  additional  question,  can  orbital 
maintenance  be  optimized  if  thrusters  are  fired  at  a  fixed 
angle  to  the  local  horizon  and  if  so,  what  is  the  angle  or, 
preferably,  range  of  angles  that  achieve  optimality? 

Ross  and  Melton  [Ref.  3]  develop  their  theory  through  the 
methods  of  optimal  control  theory.  It  is  proposed  here  that, 
while  ideally  accomplished  in  that  manner,  satisfactory  and 
enlightening  results  may  be  obtained  through  the  use  of 
relatively  unsophisticated  mathematics  and  the  aid  of  computer 
modelling  techniques.  To  exemplify  this  statement,  let  us 
examine  the  mathematical  development  in  Ross  [Ref.  3:  pp. 
1-3] .  Drag  is  given  by 

D   =  -±p(z)SzefCDvy  (1) 

where  p(r)  is  atmospheric  density,  Sref  is  the  reference 
surface  area  of  the  spacecraft  upon  which  the  aerodynamic 
forces  act,  CD  is  spacecraft's  coefficient  of  drag  and  v  is 
its  velocity.  The  equations  of  motion  can  be  written  as 

t   -  Y  (2) 
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Applying  the  principles  of  Pontryagin,  the  Hamiltonian  is 
given  by 
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The  costates  are  then  written  as 
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A  closed  form  solution  to  these  equations  does  not  exist  and 
only  after  initial  and  boundary  conditions  have  been 
determined  may  the  solution  be  obtained  through  numerical 
integration.  This  is  a  cumbersome  method  for  determining  the 
optimal  direction  of  thruster  firing.  In  this  thesis,  we  look 
into  the  possibility  of  a  constant  vectoring  scheme  that  may 
result  in  nearly  identical  performance. 


A  computer  generated  model  will  strive  to  maintain 
spacecraft  trajectory  within  ±1.5  kilometers  of  the  injection 
altitude  utilizing  a  periodic  fixed  angle  transverse  thrusting 
control  scheme  whose  direction  is  maintained  at  a  fixed  angle 
relative  to  the  local  horizon.  The  propellant  consumed  will 
then  be  compared  to  that  of  the  same  satellite  employing  a 
control  law  that  sets  thrust  equal  to  drag  at  every  point 
within  the  orbit. 


III.   SIMPLIFIED  PARAMETRIC  FORMULATION 

A  method  less  elegant  than  optimal  control  theory,  but 
nonetheless  valid,  is  that  of  parametric  variation.  The 
equations  of  motion  are  developed  for  the  satellite's  orbital 
trajectory  and  certain  parameters  varied  to  achieve  "optimal" 
control  of  orbital  variations. 

A.   DEVELOPMENT  OF  THE  EQUATIONS  OF  MOTION 

For  the  purpose  of  building  a  solid  foundation,  certain 
simplifying  assumptions  are  made.  Orbital  motion  is  assumed  to 
be  coplanar,  the  initial  spacecraft  orbit  is  assumed  to  be 
circular  and  since  the  spacecraft  is  a  non-lifting  (blunt) 
body,  drag  is  the  net  aerodynamic  force  acting  upon  it. 

The  equations  of  motion  for  this  two-body  system  can  be 
written  as 

a,  =  £—  (9) 

a£r  =  E—  (10) 

where  ar  is  the  radial  acceleration  of  the  spacecraft,  atr  is 
its  transverse  acceleration,  Fr  is  the  sum  of  the  external 
forces  in  the  radial  direction,  Ftr  is  the  sum  of  the  external 
forces  in  the  transverse  direction  and  m  is  spacecraft  mass. 


The  external  forces  acting  on  the  spacecraft  are  the 
gravitational  field,  aerodynamic  forces  and  thrust.  Figure  1 
illustrates  the  coordinate  system  and  the  net  forces. 


Transverse  axis 


Satellite  orbit 


Figure  1  Graphical  Representation  of  Coordinate  System 


Referring  to  Figure  1,  it  is  clearly  seen  that  the  components 
for  drag  are  given  by 


Dz   =  -D  sin(y) 


(11) 


Dt_  =  -D  cos(y) 


(12) 


Likewise,    thrust    is    written    as 

?x  -  ^  cos  (a)  (13) 

T*    =    Tmax   Sin(tt)  (14) 

The  angle  y     is  the  flight  path  angle,  defined  by  the 

intersection  of  the  velocity  vector  and  the  transverse  axis. 

The  angle  a  is  called  the  thrust  angle,  defined  by  the 
intersection  of  the  thrust  vector  and  the  transverse  axis.  The 
equations  of  motion,  then,  become 

f  -  62r  =  — H-  -  Bl   +  1*  (15) 

r z         m         m 


6r  +  26r  =  -Biz.   +  111  (16) 

m  in 

where  |i  is  the  Earth's  gravitational  constant. 


B.   DETERMINATION  OF  THE  CONTROL  VARIABLE 

Following  the  development  of  the  previous  chapter,  it  is 
desirable  to  maintain  the  spacecraft  within  a  radial  band  of 
a  predetermined  width.  This  may  be  accomplished  by  directly 
controlling  either  radius  or  eccentricity.  Radius  is  the 
control  variable  of  choice  for  a  number  of  reasons:  it  is 
found  directly  from  integration  of  the  equations  of  motion, 
changes  are  easily  visualized  and  radius  control  provides 
indirect  control  of  the  eccentricity.  It  is  clear  that  by 
varying  the  thrust,  control  of  satellite  radius  is  possible. 
Examination  of  the  thrust  equations  (Equations  13  and  14), 
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presents  two  methods  by  which  thrust  may  vary,  changes  of 
amplitude  or  change  in  the  direction  of  the  thrust  vector. 

The  following  chapter  presents  a  method  by  which  radial 
deviation  is  controlled  through  variation  of  the  thrust  angle 
and  then  tested  by  varying  the  thrust  magnitude.  A  computer 
model  is  developed  that  simulates  the  trajectory  of  a 
spacecraft,  graphically  demonstrating  the  effects  imposed  on 
it  through  variation  of  the  direction  of  the  thrust. 
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IV.  DEVELOPMENT  AND  TESTING  OF  THE  COMPUTER  MODEL 

A.   COMPUTER  PROGRAM  DEVELOPMENT 

As  indicated  in  the  previous  chapter,  a  computer  program 
is  developed  to  simulate  spacecraft  orbital  trajectories  and 
is  included  in  Appendix  A.  The  program  is  written  in  the 
fortran  language  and  employs  a  fourth-order  Runge-Kutta 
numerical  integration  routine  to  integrate  the  equations  of 
motion.  The  program  consists  of  six  sections,  the  main  program 
and  five  subroutines.  The  main  program  controls  input  and 
output  while  the  subroutines  provide  various  other  functions 
necessary  for  accurate  simulation  of  orbital  trajectories. 

The  first  subroutine  calculates  drag  experienced  by  the 
spacecraft.  Initially,  a  model  incorporating  constant 
atmospheric  density  is  used.  This  facilitates  verification  of 
the  program,  after  which,  an  exponential  density  model  is 
used.  It  is  acknowledged  that  more  accurate  atmospheric 
density  models  exist,  however,  the  exponential  model  provides 
satisfactory  accuracy  over  the  range  of  travel  experienced  by 
the  satellite  (±  1.5  kilometers  from  initial  orbit  R0)  as 
controlled  by  the  simulation.  The  second  subroutine  contains 
the  equations  of  motion  governing  the  spacecraft's  orbital 
trajectory.  To  facilitate  handling,  the  equations  were  broken 
into  parts.  Solving  Equations  (15)  and  (16)  for  acceleration, 
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it  is  seen  that  the  right-hand  side  of  the  radial  equation  has 
four  components  and  the  angular  equation  has  three.  These 
components  are  labeled  A,  B,  C,  and  E  for  the  radial  equation 
and  P,  Q  and  S  for  the  equation  governing  angular  motion.  The 
third  subroutine  contains  the  fourth-order  Runge-Kutta 
numerical  integration  routine  used  to  integrate  the  equations 
of  motion.  The  next  subroutine  calculates  the  parameters  of 
the  satellite's  osculating  orbit.  The  last  subroutine  contains 
the  control  law  governing  the  activation  and  deactivation  of 
the  thrusters  responsible  for  the  periodic  maintenance  of  the 
satellite's  orbit. 

The  following  parameters  define  the  specifications  around 
which  the  computer  model  was  developed. 

•  Spacecraft  mass  =  20,000  kg. 

•  Spacecraft  frontal  area  (Sref)  =  60  m2 . 

•  Coefficient  of  drag  (CD)  =  2.2. 

•  Altitude  of  spacecraft's  orbit  (Rc)  =  260  km. 

•  Atmospheric  density  at  R0  (p)  =  8.3130  x  10"11  kg/m3. 

B.   PROGRAM  VALIDATION 

Program  development  proceeded  in  stages,  with  each  stage 
requiring  validation  prior  to  beginning  the  next.  Initially, 
all  external  forces  except  gravitation  were  neglected. 
Clearly,  radius,  speed,  angular  momentum  and  specific  energy 
must  remain  constant  for  the  program  to  be  considered  valid. 
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Having  accomplished  that  elementary  stage,  phase  two 
introduced  drag  (Equation  (1)).  Taking  into  account  the  fact 
that  aerodynamic  forces  acting  on  a  spacecraft  are  very  small, 
and  consequently,  the  time  required  for  significant  changes  to 
occur,  very  large,  atmospheric  density  was  increased  by 
approximately  three  orders  of  magnitude  in  order  to  reduce 
computer  run  time. 

Initially  verifying  that  radius  continually  decreased,  the 
accuracy  was  tested  by  comparing  the  difference  in  altitude 
per  orbit  calculated  by  the  program  to  that  calculated 
manually  by  simplified  equations.  This  is  accomplished  through 
the  use  of  a  "rule-of-thumb"  and  is  developed  in  the  following 
section . 

1 .   Development  of  a  Rule  of  Thumb 

Work  done  by  drag  is  a  function  of  the  path  travelled 
by  the  spacecraft.  Therefore,  the  amount  of  work  done 
corresponds  to  the  change  in  energy  of  the  spacecraft,  which 
is  given  by 

Work  done  =  &E  =  Drag  x  2nr  (17) 

Spacecraft  specific  energy  is  given  by 

r=  *L  -  J£  =  — £-  (is) 

2    r  2a 

Total  energy  E  is  equal  to  the  specific  energy  multiplied  by 
spacecraft  mass.  Assuming  that  the  spacecraft  is  in  a  circular 
orbit,   the  semi-major  axis  a  is  equal  to  the  radius  r. 
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Performing  this  substitution  then  differentiating  both  sides 
of  the  equation  while  holding  mass  constant  yields 

dE  =  -ML  dr  (19) 

2r2 

Setting  Equations  (17)  and  (19)  equal  to  each  other  and 
solving  for  the  change  in  radius  yields 

Ar  =  47rr3'P  (20) 

Equation  (20)  is  the  decrease  in  radius  per  orbit  of  the 
spacecraft  due  to  drag.  Despite  inaccuracies  resulting  from 
the  simplifying  assumptions,  this  rule  of  thumb  is  accurate  to 
within  a  few  percent.  The  difference  in  spacecraft  radius 
calculated  by  the  computer  program  matched  that  of  the  thumb 
rule  within  a  few  percent  thereby  validating  the  model  through 
this  point  in  its  development.  Phase  three  introduced  thrust 
while  setting  drag  equal  to  zero.  Clearly,  any  results  other 
than  steadily  increasing  radius,  angular  momentum,  and 
specific  energy  would  have  been  cause  for  program 
invalidation . 

2 .   Development  of  a  Control  Law 

The  purpose  of  the  control  law  is  to  maintain 
satellite  radius  within  a  prespecified  bandwidth.  By 
monitoring  certain  variables,  activation  and  deactivation  of 
station  keeping  thrusters  can  be  determined.  Keeping  the 
control  law  simple  in  order  that  computer  memory  and  run  time 
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related  to  station  keeping  be  kept  to  a  minimum  is  also  a 
desirable  goal.  With  these  facts,  control  law  design  proceeded 
as  described  in  the  next  section. 
3 .   Control  Law  design 

a.  Approximate  method 

Specific  energy  is  a  function  of  spacecraft  radius. 
For  that  reason,  specific  energy  is  a  very  useful  parameter 
that  can  be  used  to  maintain  the  satellite  within  the 
specified  bandwidth.  The  energy  lost  when  the  satellite's 
radius  decreases,  is  dependent  upon  the  path  taken  by  the 
satellite  in  its  descent  as  illustrated  in  Figure  2  below 
(exaggerated  for  clarity)  . 


Figure  2  Path  Travelled  By  A  Spacecraft  Between  Two  Orbits 


If  S  represents  the  path  travelled  by  the  spacecraft,  then  Sx 
is  its  initial  position  and  S2  its  final  position.  The  force 
acting  against  spacecraft  motion  is  drag,  which  directly 
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opposes  the  velocity  vector.  With  this  in  mind,  the  energy 
loss  is  given  by 

LE  =  [D  dS  (21) 

The  arc  length  S  is  a  function  of  radius  and  angle  turned 
through,  and  is  given  by 

S  =  r6  (22) 

Differentiating  the  above  equation  yields 

dS  =  rdd   +  6dr  (23) 

Substituting  Equation  (23)  into  Equation  (21)  yields  energy 
loss  in  terms  of  the  known  variables,  r  and  9 


AE  =  fnddr   +  JDrdd  (24) 


*0 

Employing  Simpson's  rule,  the  loss  in  energy  can  be 
approximated  fairly  closely.  Assuming  atmospheric  density, 
flight  path  angle,  thrust  angle  and  thrust  are  all  constant, 
Simpson's  rule  applied  over  ten  iterations  yields  the  energy 
lost  in  moving  from  Rc  to  Rf 


AB  = 


Rf M^eeo)  +  -^-JDr(29)  <25> 

^   30  )  30 


Referring  to  Figure  2,  Rf  is  the  point  at  which  the  thrusters 
will  fire.  In  terms  of  the  control  law,  Rf  is  the  control 
variable  and  can  be  relabelled  as  AR  since  it  is  variable  and 
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determines  the  width  of  band  in  which  the  satellite  is 
maintained.  As  developed  for  this  model,  AR  is  the  initial 
radius  R0  minus  one  kilometer.  This  value  is  chosen  so  as  to 
maintain  maximum  orbital  deviation  within  ±  1.5  kilometers  of 
R0.  Now  that  the  energy  loss  has  been  approximated,  the 
objective  is  to  determine  the  length  of  time  to  fire  the 
thrusters  in  order  to  replace  the  energy.  The  return  path  of 
the  spacecraft  is  a  function  of  the  thrust. 


A£  =  f  (f  -  D)    •  ds  (26) 


Recognizing  that  ds_  is  related  to  the  time  rate  of  change  of 
the  position  vector  or  arc  length  S_,  Equation  (26)  becomes 


LE  -    f(f  -  D)»v  dT 

Reducing   Equation    (27)    to    component    form   yields 

AE  =  f(Tz-Dr)vrdt  +  |(Ttr  -  DtI)vtIdt 


(27) 


(28) 


Integrating  and  solving  for  At  yields  the  length  of  time  that 
the  thrusters  must  fire  in  order  to  replace  the  energy  lost 
due  to  drag. 

At  =  — (29) 

(rr  -Dz)vz  +    [Ttz  ~  DtI)vtI 

Examination  of  Figure  1  reveals  that  velocity  can  be  written 
as 
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v  =  vsinCy)  (30) 

and 

vtI   =  vcos  (y)  (31) 

Figures  1  through  6  of  Appendix  B  illustrate 
results  obtained  from  this  method  for  thrust  angles  of  60°  and 
65°.  Examination  of  plots  of  thrust  versus  orbit  (Figures  2 
and  5)  reveal  the  inadequacy  of  this  model.  Thruster  firing 
times  are  seen  to  be  on  the  order  of  orbits  rather  than 
fractions  of  orbits.  This  is  due  to  two  factors,  first  and 
most  important  being  that  the  energy  change  is  calculated  from 
the  initial  orbit  to  the  point  AR  where  the  thrusters  begin 
firing.  The  problem  arises  from  the  fact  that  the  satellite  is 
still  in  a  descent  at  this  point  and  continues  to  descend 
until  its  motion  is  reversed  through  the  opposing  force  of  the 
thrusters.  As  a  result,  the  satellite  loses  more  energy  than 
is  replaced.  The  second  problem  arises  from  the  inaccuracies 
inherent  in  the  assumptions  required  to  perform  the 
approximation.  While  the  first  problem  is  easily  resolved,  the 
changes  in  atmospheric  density  and  flight  path  angle,  while 
very  small,  are  not  constant  and  the  resulting  errors  combine 
to  render  this  model  unsatisfactory.  A  more  accurate  method  is 
to  calculate  the  energy  loss  directly  using  the  variables 
derived  from  integrating  the  equations  of  motion. 
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b.    Direct  method 

Using  Equation  (18) ,  the  initial  and  instantaneous 
energies  may  be  calculated  at  any  instant  during  the 
trajectory  of  the  spacecraft.  As  before,  thrusters  will  fire 
when  the  spacecraft  descends  below  AR.  The  program  then 
calculates  the  spacecraft's  specific  energy  each  time  the 
equations  of  motion  are  integrated.  Comparing  this  value  to 
that  of  the  initial  circular  orbit,  the  control  law  commands 
the  thrusters  to  fire  until  they  are  equal.  Examination  of 
Figures  1  through  6  of  Appendix  C  reveals  that  while  results 
are  closer  to  those  expected,  this  method  also  appears  to  be 
inadequate.  Thruster  firing  times  (Figures  2  and  5)  remain 
excessively  long.  Reevaluation  of  the  model  suggested  that  the 
solution  might  be  a  function  of  the  thrust  to  drag  ratio. 
Increasing  thrust  to  a  value  of  300  N  and  then  examining 
results  for  thrust  angles  of  60°,  65°  and  70°  revealed  this  to 
be  the  solution.  Results  are  found  in  Appendix  D.  Thrust 
angles  of  65°  and  70°  maintain  radial  deviation  within  the 
prescribed  bandwidth  of  three  kilometers  with  a  trend  that 
indicates  they  will  remain  so  indefinitely.  Thrust  plots 
(Figures  2,  5  and  8)  illustrate  that  the  thrusters  are  firing 
over  a  small  portion  of  an  orbit  instead  over  a  period  of  many 
orbits  as  before.  Plots  of  spacecraft  radius  versus  orbit  are 
included  in  each  appendix  to  illustrate  the  success  of  the 
control  law  in  maintaining  radial  deviation  within  the  three 
kilometer   band.   Energy   plots   are   also   included   as 
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corroberating  illustrations  of  the  spacecraft's  energy  level 
at  each  point  in  its  trajectory. 
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V.  ANALYSIS  AND  RESULTS 

Analysis  of  the  validated  model  is  performed  in  four 
steps.  Prior  to  testing,  drag  was  returned  to  its  true  value 
of  8.3130  x  10"::  kg/m3  and  thrust  was  set  equal  to  a  value  of 
25  N.  This  value  is  determined  by  multiplying  the  thrust  to 
drag  value  of  the  test  model  by  the  drag  actually  experienced 
by  the  satellite.  The  thrust  to  drag  ratio  with  the  increased 
density  is 

300N      =75.834  (32) 


3  .956N 


With  p  =  8.313  x  10"::  kg/m3,  drag  is  calculated  to  be  0.329  N. 
Multiplying  this  to  the  value  in  Equation  (32)  yields  a  thrust 
of  25  N. 

A.   USE  OF  A  CONSTANT  DENSITY  MODEL 

Initial  testing  maintained  a  constant  density  while 
varying  thrust  angle.  The  results  obtained  from  these  tests 
are  interesting.  As  illustrated  in  Appendix  E,  approximately 
six  and  a  half  orbits  are  required  for  the  spacecraft's  radius 
to  decay  to  AR.  Thruster  firing  places  the  spacecraft  into  a 
slightly  elliptical  orbit  (typically,  eccentricities  were 
found  to  be  on  the  order  of  10"4)  .  Examination  of  energy  plots 
reveal  that  the  energy  of  the  elliptical  orbit  is  very  close 
to  that  of  the  initial  circular  orbit.  This  results  in  one 
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large  energy  change  as  the  spacecraft  initially  descends  from 
Rc  and  is  returned,  then,  subsequent  small  changes  are 
required  once  the  spacecraft  is  established  in  its 
"elliptical"  orbit.  A  thrust  angle  of  70°  is  seen  to  yield  the 
best  results  (see  Figure  6  of  Appendix  E)  .  Radial  deviation 
remains  just  within  the  specified  radial  band  with  a  trend 
that  indicates  it  will  do  so  indefinitely.  At  angles  of  less 
than  70°,  radial  deviation  steadily  increases  until  it  exceeds 
the  prescribed  limits.  An  illustration  is  provided  in  Appendix 
E,  for  a  thrust  angle  of  65°.  Thrust  angles  greater  than  70° 
follow  a  trend  illustrated  by  Figure  7  of  Appendix  E  (01=75°) 
until  reaching  approximately  85°.  Above  85°,  increasingly 
larger  values  of  thrust  are  required  to  maintain  the 
spacecraft  within  the  prespecified  radial  band. 

B.   USE  OF  AN  EXPONENTIAL  ATMOSPHERIC  MODEL 

Use  of  a  constant  atmospheric  density  model  permitted  the 
determination  of  an  optimum  angle  that  could  be  compared  with 
that  determined  by  a  more  accurate  exponential  atmospheric 
density  model.  Also,  realizing  that  density  changes  would  be 
small  within  a  three  kilometer  band,  the  constant  density 
model  divulged  a  reasonable  facsimile  of  results  obtained  from 
the  exponential  model.  As  stated  previously,  more  accurate 
density  models  exist,  for  example,  the  Jacchia  Atmospheric 
Model   (J70) .   The   exponential   model,   however,   provides 
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sufficient  accuracy  for  the  development  and  testing  done  in 
this  thesis.  With  this  in  mind,  density  is  now  given  by 


P  =  Poe 


-P(r-r0)  (33) 


where  (3  is  determined  to  be  -2.12  x  10"5  irf1  [Ref.  7]. 
As  predicted,  the  results  are  nearly  identical  to  those 
obtained  from  the  constant  density  model.  A  thrust  angle  of 
70°  maintains  the  spacecraft  within  the  three  kilometer  band 
with  the  same  trends  as  described  in  the  previous  section. 
Results  are  illustrated  in  Appendix  F. 

1 .   Comparison  with  a  Forced-Keplarian  Model 

It  is  now  possible  to  compare  results  obtained  from 
this  model  with  that  of  a  spacecraft  experiencing  forced- 
Keplarian  motion  (thrust  equal  to  drag,  resulting  in  a 
circular  motion) .  Modifications  to  the  program  to  obtain  a 
model  in  which  thrust  is  equal  to  drag  are  very  simple.  The 
drag  subroutine  calculates  drag  then  sets  thrust  equal  to  it. 
The  subroutine  containing  the  thrust  control  law  is  removed 
from  the  program  entirely  since  thrusters  will  fire 
continuously  as  long  as  propellant  is  available.  Drag  always 
opposes  the  velocity  vector  which,  as  previously  shown,  is 
defined  by  the  flight  path  angle  y.  If  the  thrust  is  exactly 
opposite  and  equal  to  the  drag  force,  then  the  thrust  angle  a 
is  equal  to  the  flight  path  angle.  Therefore,  modifications 
consisted  of  setting  a  =  y,    removing  the  subroutine  containing 
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the  thrust  control  law  and  setting  thrust  equal  to  the  drag 
calculated  in  the  appropriate  subroutine.  Since  the  spacecraft 
is  initially  in  a  circular  orbit,  validation  of  these 
modifications  is  achieved  when  the  spacecraft's  orbital 
parameters  remain  unchanged  over  the  test  period;  in  this 
case,  20  orbits. 

Finding  that  the  model  performed  as  expected,  plots  of 
propellent  consumed  over  the  test  period  are  compared  with 
those  for  the  spacecraft  experiencing  non-Keplarian  motion  at 
a  thrust  angle  of  70°.  Results  are  contained  in  Appendix  G.  It 
is  clearly  seen  that  orbital  maintenance  using  forced- 
Keplarian  motion  is  superior  to  that  using  non-Keplarian 
motion . 

2 .   Further  Testing 

To  determine  the  "robustness"  of  the  model,  two 
additional  tests  were  performed.  In  the  first,  thrust  was 
maintained  at  25N  while  specific  impulse  was  varied  over  the 
range  of  200  sec,  valid  for  chemical  reaction  engines,  to  2000 
sec  which  is  valid  for  electrothermal  engines.  In  the  second 
test,  specific  impulse  is  returned  to  its  initial  value  of  300 
sec  while  thrust  is  varied  over  a  range  extending  from  IN  to 
35N. 

a.    Constant   Thrust,    Varied  Specific  Impulse 

Varying  specific  impulse  while  holding  thrust 
constant,  forces  the  rate  of  fuel  consumption  to  change. 
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Specific  impulse  is  given  by  given  by 

ISP   =  -3T-  <34> 

where  I  is  the  specific  impulse,  T  is  the  magnitude  of  the 
thrust,  g0  is  the  gravitational  acceleration  at  sea-level  and 
m  is  the  change  in  mass. 

Clearly,  decreasing  thrust  results  in  a  corresponding 
decrease  in  the  rate  of  propellant  consumed  over  a  given 
period  of  time.  This  is  graphically  represented  in  Appendix  G. 
The  absolute  quantity  of  fuel  consumed  decreases  with 
increasing  specific  impulse,  but  the  percent  difference 
between  the  reboost  and  forced  Keplarian  methods  remain 
virtually  unchanged." 

The  percent  difference  is  calculated  by  taking  the  ratio 
of  the  values  of  propellant  mass  at  a  specific  time  for 
thrust-equal-drag  plots  and  reboost  plots  to  determine  the 
relative  change  between  them.  This  provides  a  truer  comparison 
of  the  two  schemes  than  does  simply  comparing  the  end  values 
of  the  plots. 

The  results  illustrated  in  Appendix  G  indicate  that,  while 
playing  a  significant  role  in  propulsion  system  optimization, 
specific  impulse  is  not  a  function  of  the  method  used  to 
perform  orbital  maintenance  and  will  not  affect  the  particular 
outcome  of  one  method  more  than  another. 
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b.    Constant   Specific  Impulse,    Varied  Thrust 

The  previous  section  illustrated  the  significance 
played  by  specific  impulse  in  minimizing  fuel  expenditure 
during  orbital  maneuvers.  Obviously,  engine  "size"  plays  an 
equally  important  role  in  that  process.  It  is  expected  that  as 
thruster  size  increases,  fuel  expenditure  will  increase. 
Appendix  H  bears  this  out.  While  it  is  theoretically  a  simple 
matter  to  choose  the  proper  specific  impulse  (bigger  is 
better),  this  is  not  the  case  when  choosing  thruster  size. 
Examination  of  the  plots  in  Appendix  H  reveal  that  when 
comparing  engines  over  a  certain  time  span,  the  end  results  do 
not  provide  a  ready  solution.  Figure  1  is  a  case  in  point. 
Although  this  case  (IN  thruster,  Isp=300  sec)  results  in  the 
least  amount  of  fuel  expended  over  the  range  of  thrusts 
chosen,  it  is  obvious  that  this  is  not  a  wise  choice  of 
thruster  size.  The  thruster  fires  continuously  from  its 
initiation  until  the  end  of  the  test  period.  The  thruster  is 
clearly  too  small.  A  five  newton  thruster,  i.e.,  Figure  2, 
seems  to  be  a  viable  engine  size,  although,  without  more 
information,  it  is  difficult  to  determine  positively. 

Relatively  large  thrusters  burn  for  shorter  periods 
of  time  than  smaller  thrusters  to  achieve  a  common  result,  but 
each  burn  expends  more  fuel.  Smaller  thrusters,  on  the  other 
hand,  expend  less  fuel  for  a  given  burn  time,  yet  must  burn 
longer  to  achieve  the  same  results.  The  fundamental  result  of 
this  test  is  that  thruster  size  is  a  tradeoff  variable  that  is 
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to  be  used  in  conjunction  with  other  factors  to  achieve 
desired  results,  such  as  minimization  of  thruster  burn  time 
during  orbital  maintenance. 
3.   Final  Tests 

Previous  sections  have  shown  that  within  a  narrow 
radial  bandwidth,  a  model  using  thrust  equal  to  drag  is 
superior  to  one  using  fixed  thrust-angle  reboost  techniques. 
The  question  arises  as  to  whether  these  results  will  remain 
valid  for  larger  bandwidths. 

As  described  earlier,  the  control  law  commands  the 
thrusters  to  fire  when  the  spacecraft  orbit  has  decayed  a 
certain  distance  below  the  reference  orbit.  Results  are 
examined  for  cases  where  the  spacecraft  is  allowed  to  descend 
20  km,  30  km  and  40  km  below  R0.  Thrust  is  fixed  at  300  N 
while  specific  impulses  vary  between  300  sec,  500  sec  and  2000 
sec.  As  for  the  case  of  a  3  km  bandwidth,  a  thrust  angle  of 
70°  maintained  the  spacecraft  within  the  desired  bandwidth  and 
was  therefore  used  for  all  cases  described  below.  The  tests 
are  performed  over  a  period  of  100  orbits.  In  order  to  reduce 
computer  time,  atmospheric  density  was  once  again  increased  to 
a  value  of  1  x  10-9  kg/m3. 

a.  Case  1:    Thrusters  Fire  20  km  Below  R0 

Rather  than  choosing  a  specific  bandwidth  and 
adjusting  the  control  law  to  achieve  it,  the  spacecraft's 
orbit  was  allowed  to  decay  a  certain  distance  prior  to 
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activation  of  the  control  thrusters  and  the  resulting 
bandwidth  measured.  This  provided  expediency  since  the  control 
law  determined  the  bandwidth  rather  than  having  to  be  adjusted 
to  achieve  it.  The  results  are  the  same  in  either  case  so  no 
accuracy  is  lost  with  this  method.  The  figures  in  Appendix  I 
illustrate  the  results  of  this  case. 

Allowing  the  spacecraft's  orbit  to  decay  20  km 
prior  to  activation  of  the  control  jets  provided  a  radial  band 
of  55  km.  Plots  of  expended  propellant  mass  versus  orbit  are 
provided  for  the  three  test  cases  (Isp  =  300  sec,  500  sec  and 
2000  sec)  .  As  in  previous  cases,  thrust  equal  to  drag  provides 
a  straight  line  while  reboost  is  actually  a  series  of  steps. 
Each  step  is  a  cycle  wherein  the  spacecraft  descends  20  km  at 
which  time  reboost  occurs  (vertical  portion  of  plot) ,  after 
deactivation  of  the  thrusters,  the  spacecraft  once  again 
descends  to  the  point  where  reboost  reoccurs.  This  is 
indicated  by  the  horizontal  portion  of  the  plot  since  no  fuel 
is  being  expended  during  this  portion  of  the  trajectory.  As 
before,  the  reboost  maneuver  puts  the  spacecraft  into  a 
slightly  eccentric  orbit  (on  the  order  of  10"5)  .  This  accounts 
for  the  periodic  motion  and  high  number  of  thruster  firings 
indicated  by  the  mass  plots  in  Figures  2  through  4.  As 
expected,  increasing  specific  impulse  reduces  the  amount  of 
propellant  expended  over  the  test  cycle. 
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b.    Case   2:    Thrusters  Fire   30  km  Below  R0 

Requiring  the  control  law  to  fire  30  km  below  the 
reference  orbit  provides  a  radial  band  of  78  km.  Figure  5  is 
an  illustration  of  spacecraft  radius  over  time.  Because  AR  is 
so  large,  the  reboost  maneuver  must  occur  twice  before  the 
spacecraft  settles  into  a  periodic  trajectory  that  carries  it 
the  full  width  of  the  radial  band.  Eccentricity  of  the 
osculating  orbit,  however,  remains  on  the  order  of  10"4.  Once 
in  its  periodic  trajectory,  results  are  very  similar  to  those 
of  case  1.  Figures  6  through  8  illustrate  propellant  mass 
consumption  over  the  test  period.  The  first  two  thruster 
firings  are  clearly  evident.  Once  in  its  periodic  trajectory, 
the  mass  plots  are  very  similar  to  those  of  case  1  and  occur 
for  the  same  reasons. 

Comparisons  of  the  mass  plots  of  case  2  to  case  1 
reveal  interesting  results.  Percent  difference  comparisons  of 
case  2  to  case  1,  for  identical  specific  impulses,  provides  an 
indication  of  established  trends  from  which  inferences  of 
future  results  might  be  made. 

We  see  that  for  case  1,  for  each  variation  of 
specific  impulse,  reboosting  the  spacecraft  requires  482 
percent  more  fuel  than  the  use  of  thrust  equal  to  drag 
techniques.  Case  2  reboost  maneuvers  require  578  percent  more 
fuel  than  thrust-cancel-drag  maneuvers. 
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c.    Case   3:    Thrusters  Fire   40  km  Below  R0 

Allowing  the  spacecraft's  orbit  to  decay  40  km 
prior  to  initiating  reboost  sequences  provides  a  radial  band 
of  100  km  in  which  the  satellite  trajectory  is  maintained. 
Figure  9  illustrates  radial  deviation  of  the  spacecraft.  The 
increase  in  AR  coupled  with  a  fixed  thrust  requires  the 
satellite  to  perform  three  reboost  maneuvers  before  the 
familiar  periodic  trajectory  is  achieved.  Thruster  firing  is 
clearly  evident  in  the  first  two  incidences,  as  is  the  ensuing 
orbital  decay  of  the  resulting  (slightly)  eccentric 
trajectories . 

Calculating  the  mass  expenditure  percentages 
reveals  that  reboosting  the  satellite  requires  481  percent 
more  fuel  than  does  setting  thrust  equal  to  drag.  Although 
thrust  equal  to  drag  is  still  superior,  the  difference  between 
the  two  is  decreasing.  To  further  test  this  result, 
percentages  were  calculated  for  points  43  orbits  and  97  orbits 
into  the  test  period.  All  values  were  less  than  corresponding 
values  calculated  in  case  2.  While  these  results  do  not 
provide  conclusive  evidence,  we  may  conjecture  that  a  trend  is 
developing,  indicating  that  at  some  point  results  from 
periodic  orbital  maintenance  will  equal  those  from  forced 
Keplarian  motion,  or  as  the  theory  predicts,  the  plots  will 
reverse  and  periodic  thrusting  will  become  superior. 


31 


C.   SUMMARY  OF  RESULTS 

Initial  testing,  utilizing  a  constant  atmospheric  density 
model  provided  a  baseline  against  which,  further  testing  could 
be  compared.  A  thrust  angle  of  70°  was  found  to  produce  the 
desired  results.  Fixing  the  thrust  vector  at  this  angle 
maintained  spacecraft  orbital  deviation  within  a  three 
kilometer  band  nearly  indefinitely.  Upon  determining  this 
angle,  the  computer  program  discarded  the  constant  density 
model  and  incorporated  an  exponential  atmospheric  density 
model . 

This  simulation  was  then  compared  to  a  model  in  which 
thrust  canceled  drag.  According  to  the  optimal  control  theory 
developed  in  Chapter  III,  the  simulation  (reboost  model) 
should  prove  superior  to  a  thrust-cancel-drag  model  (relative 
to  the  problem  of  fuel-minimization) .  In  fact,  reboost 
required  significantly  more  propellant  to  maintain  the 
satellite  orbit  within  the  three  kilometer  band  than  did 
thrust  equal  to  drag.  To  test  the  robustness  of  these  result, 
specific  impulse  was  varied  between  200  and  2000  sec  and 
thrust  was  varied  between  1  and'  35  N.  The  results  remained 
unchanged.  To  further  test  the  results,  the  simulation  was 
tested  over  much  wider  radial  bands.  The  results  still  proved 
thrust-cancel-drag  trajectories  superior  to  reboost  models 
although  the  difference  in  efficiency  seemed  to  decrease. 


32 


VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

As  stated  previously,  optimal  control  theory  states  that 
orbital  maintenance  using  a  technique  in  which  thrust  cancels 
drag,  is  not  optimal.  This  means  that  some  scheme  using 
periodic  thrusting  must  then  be  optimal.  Through  the 
complicated  techniques  of  optimal  control  theory,  a  thrust 
vectoring  scheme  is  shown  to  indeed  be  the  optimum.  If  the 
thrust  vector  always  points  along  the  primer  vector,  the 
trajectory  is  optimal.  A  sub-optimal  scheme  using  fixed-angle 
thrusting  and  parametric  variation  is  presented  here  as  a 
simplified  method  of  determining  the  optimality  of  orbital 
maintenance. 

In  each  series  of  tests,  minimization  of  propellent  mass 
using  fixed-angle  thrusting  has  proven  to  be  inferior  to  that 
in  which  thrust  is  set  equal  to  drag.  At  first  glance  these 
results  appear  to  contradict  the  theory  developed  by  Ross  and 
Melton  [Ref.  3].  For  small  perturbations  f orced-Keplarian 
motion  proved  to  be  superior  to  periodic  fixed-angle 
thrusting.  As  the  perturbations  increased  (indicated  by  the 
increasing  size  of  the  radial  band) ,  it  would  seem  reasonable 
to  expect  that  the  difference  in  fuel  consumption  between 
these  two  techniques  would  increase.  However,  the  results  of 
tests  described  in  section  A. 4  of  Chapter  V  reveal  that  for 
large  radial  bands,  the  percent  difference  in  propellant  mass 
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expended  between  methods  of  orbital  maintenance  using  periodic 
thrusting  and  forced  Keplarian  motion,  actually  appears  to 
decrease.  Based  on  these  results,  we  may  conjecture  that  the 
percent  difference  between  the  two  methods  tested  here  will 
continue  to  decrease  until  periodic  thrusting  yields  results 
superior  to  those  for  forced  Keplarian  motion.  Further  testing 
is  required  before  the  analytical  theory  proposed  in  [Ref .  3] 
may  be  conclusively  verified. 

The  problem  as  presented  here  is  that  of  fuel-minimization 
during  orbital  maintenance.  Solving  the  Mayer  optimality 
problem  derived  in  Ross  [Ref.  3]  yields  the  primer  vector. 
This  is  a  very  cumbersome  method  requiring  solution  of  a  two- 
point  boundary  value  problem.  If  periodic  thrusting  is  done 
along  the  primer  vector,  fuel  will  be  optimized.  This  thesis 
has  proposed  a  simpler  method  using  the  energy  balance  of  the 
satellite  to  achieve  similar  results.  Results,  however, 
indicate  that  for  small  perturbations,  forced-Keplarian  motion 
will  provide  the  best  results. 

Ultimately,  the  goal  of  this  thesis  is  to  provide  a  method 
of  fuel-minimization  that  is  practical  and  may  be  applied  to 
existing  systems.  Propulsion  systems  utilizing  vectored  thrust 
are  highly  complex  and  have  a  correspondingly  higher  chance  of 
failure.  Additionally,  the  extreme  size  of  the  perturbations 
required  before  periodic  orbital  maintenance  would  become  more 
economical  than  forced-Keplarian  motion  is  impractical. 
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Based  on  these  conclusions,  and  with  the  added  knowledge 
that  a  variable-thrust  propulsion  system  capable  of  operating 
continuously  over  the  lifetime  of  a  satellite  may  also  be 
impractical,  it  is  recommended  that  further  testing  of 
periodic  fixed-angle  transverse  thrusting  schemes  for  small 
perturbations  be  accomplished.  It  is  recommended  that  the 
Mayer  optimality  problem  described  in  Ross  [Ref.  3]  be  solved 
and  the  primer  vector  found.  The  results  should  then  be 
compared  to  those  described  in  this  thesis  to  determine  the 
actual  amount  of  savings  acquired  through  optimization.  It  is 
possible  that  the  amount  of  fuel  saved  may  not  warrant  the 
cost  of  building  a  propulsion  system  capable  of  periodic 
primer  thrusting.  Additionally,  a  comparison  of  orbital 
maintenance  techniques  presently  in  use,  with  results  found  in 
this  thesis,  should  be  accomplished  to  determine  their 
relative  efficiency. 
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APPENDIX  A 


c 
c 

C      OBJECTIVE:  DETERMINATION  OF  FIXED  THRUST  ANGLE  TO  MAINTAIN 

C  ORBITAL  DEVIATION  WITHIN  A  PREDETERMINED 

C  BANDWIDTH . 

C 

C 

C      VARIABLE  DEFINITIONS: 

C 

C  X(l)  =  RADIUS  (METERS) 

C  X(2)  =  RADIAL  VELOCITY  (METERS  PER  SECOND) 

C  X(3)  =  THETA   (RADIANS) 

C  X(4)  =  ANGULAR  VELOCITY  (RADS  PER  SECOND) 

C  XDOT(l)  =  TIME  DERIVATIVE  OF  X ( 1 ) 

C  XDOT(2)  =  TIME  DERIVATIVE  OF  X(2) 

C  XDOT(3)  =  TIME  DERIVATIVE  OF  X(3) 

C  XDOT(4)  =  TIME  DERIVATIVE  OF  X(4) 

C  RO  =  REFERENCE  ORBIT 

C  D  =  DRAG  (N) 

C  EO  =  SPECIFIC  ENERGY  OF  REFERENCE  ORBIT 

C  E  =  SPECIFIC  ENERGY 

C  MU  =  GRAVITATIONAL  CONSTANT 

C  M  =  SPACECRAFT  MASS  (KG) 

C  GAMMAR  =  FLIGHT  PATH  ANGLE  (RADIANS) 

C  GAMMAD  =  FLIGHT  PATH  ANGLE  (DEGREES) 

C  TH  =  THRUST  (N) 

C  TMAX  =  BLOWDOWN  (MAXIMUM)  THRUST 

C  ALPHAR  =  THRUST  ANGLE  (RADIANS) 

C  CD  =  COEFFICIENT  OF  DRAG 

C  RHOO  =  REFERENCE  ATMOSPHERIC  DENSITY 

C  RHO  =  CALCULATED  ATMOSPHERIC  DENSITY 

C  SPI  =  SPECIFIC  IMPULSE 

C  V  =  VELOCITY 

C  SREF  =  REFERENCE  SURFACE  AREA 

C  GO  =  GRAVITATIONAL  ACCELERATION 

C  H  =  INCREMENT  OF  TIME  (STEP  SIZE) 

C  PTI  =  PRINT  TIME  INTERVAL  (STEP  SIZE) 

C  T  =  BEGIN  TIME 

C  TF  =  FINAL  (END)  TIME 

C  a  =  SEMI -MAJOR  AXIS 
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C  e    =    ECCENTRICITY 

C 

c 

C      START  PROGRAM 

PROGRAM  ACTORB 
C 
C 
C      VARIABLE  DECLARATION 

IMPLICIT  REAL*8 (A-H, L-Z) 

DIMENSION  X(4),XDOT(4) 
C 
C 
C      DEFINITION  OF  CONSTANTS 

PI=3.14159265359D+0 

G0=9.806D+0 

Tl=0D+0 

N=4 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C      MAIN  PROGRAM 

c 

OPEN(10,FILE=' INIT' , STATUS='OLD' ) 

OPEN(ll,FILE  =  'OUT'  , STATUS='NEW  ) 

OPEN(12,FILE='ORBEL' , STATUS='NEW ) 

OPEN(13,FILE='RAT'  ,STATUS  =  'NEW  ) 
C 

READ (10, 1)R0, V0,M,TMAX,T,TF,H,PTI,CD,MU,RHO0,SREF,SPI 
1      FORMAT (13 (/,21X,D13.7) ) 
C 
C 

PRINT*, 'ENTER  ALPHA' 

READ*,DEG 

ALPHAR=DEG*P I / 1 8  0 . 0D  +  0 
C 
C 

INDEX=0 

KOUNT=l 

X(1)=R0 

X(2)=0D+0 

X(3)=0D+0 

X(4)=1.1673449D-3 

E0=M* ( (V0*V0) /2-MU/R0) 
C 
C 

WRITE (11, *), '    TIME   ORBITS    RADIUS   VELOCITY   ALPHA 
*ANGLM  ENERGY       TMAX' 

WRITE (11, *), '    (sec)  (km)    (km/sec)   (deg) 

*' 

WRITE (12, *), '       T         a         e      APOGEE     PERIGEE 
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♦PERIOD  ' 

WRITE (13, *),  '    TIME    ORBITS       DRAG        TH 
*  MASS         GAMMA' 

C 
C 

C      CALCULATIONS 
C 
C 
100    CALL  DRAG(SREF,CD,X,R0,RHO0,D,T, V) 

CALL  ONOFF (R0,X,TMAX, SPI , GO , TH, E0, GAMMAR, ALPHAR, V, D] 

CALL  DIFFEQ(X,XDOT,MU,D,M, TMAX, ALPHAR,  T,  TH) 

CALL  RK4  (T,X,XDOT,N,H,  INDEX,  T) 
C 
C 

IF  (INDEX  .NE.  0)  GO  TO  100 
C 
C 

C      UNIT  CONVERSIONS 
C 

R=X(1) /1000 

SPEED=(  ((X(2)*X(2))  +  (X(1)*X(4))**2)**0.5)/1000 

V=SPEED*1000 

M=M- ( (TH*H) / (SPI*G0) ) 

MASS=20000-M 

GAMMAR=ATAN(X(2)  /  (X ( 1 )  *X ( 4 )  )  ) 

GAMMAD=GAMMAR*180 . 0D+0/PI 

ANGM= (X ( 1 ) *V*COS (GAMMAR) ) 

ALPHA=ALPHAR*180 .0D+0/PI 

ENERGY=( (V*V) /2)-(MU/X(l) ) 

ORBITS=T/5382.458 


CALL  ORBDAT (ENERGY, ANGM, PI, a, e, APOGEE, PERIGE, PERIOD) 


C 

C 

C      OUTPUT 

C 

IF  (KOUNT   .LT.  DNINT (PTI/H) )  GO  TO  200 
C 

WRITE (11,2) T, ORBITS, R, SPEED, ALPHA, ANGM, ENERGY, TH 

WRITE (12,  3) T, a, e, APOGEE, PERIGE, PERIOD 

WRITE (13,4) T, ORBITS, D,TH, MASS, GAMMAD 
C 

2  FORMAT (2X,F7.0, IX, F5 . 2, 3X, F8 . 3, 2X, F6 . 4, IX, F6 . 1,  2X, 
*F14.2,2X,F12.2,2X,F5.0) 

3  FORMAT (2X,F7.0, IX, F10 . 3, 2X, F4 . 3, IX, F10 . 3, IX, F10 . 3, 2X, F10 . 2 ) 

4  FORMAT(2X,F7.0,2X,F5.2,2X,F12.9, IX, F5 . 0, IX, F15 . 3, IX,  F10 . 4) 
C 

C 
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KOUNT=0 
200    KOUNT=KOUNT+l 

IF  (TF  .GE.  T)  GO  TO  100 
C 

c 

END 
C 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE  DRAG (SREF, CD, X, R0, RHO0, D, T, V) 

IMPLICIT  REAL*8 (A-H,L-Z) 

DIMENSION  X(4),XDOT(4) 
C 
C 

V=(  (X(2)  *X(2)  )  +  (X(l)  *X(4)  )  **2)  **0.5 

BETA=2 . 12D-05 

RHO=RHO0*EXP (-BETA* (X (1) -R0) ) 

D=0 . 5D0*RHO*CD*SREF*V*V 
C 
C 

RETURN 

END 
C 
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

SUBROUTINE  DIFFEQ (X, XDOT,MU, D,M, TMAX, ALPHAR, T,  TH) 
IMPLICIT  REAL*8 (A-H, L-Z) 
DIMENSION  X(4),XDOT(4) 


A=X(1)  *X(4)  *X(4) 
B=MU/ (X(l) *X(1) ) 
C=(D/M)  *  (X(2)  /  (  (  (X(2)  *X(2)  )  +  (X(l)  *X(4)  )  **2)  **0.5)  ) 

E= (TH/M) *SIN (ALPHAR) 


P=2*X(4)  *X(2)  /X(l) 

Q=(D/(X(1)*M)  )*(  (X(1)*X(4)  )/(((X(2)*X(2)  )  +  (X(l)*X(4)  )**2)**0.5)) 

S=(TH/ (X(l) *M) ) *COS (ALPHAR) 


XDOT(l)=X(2) 

XDOT(2)=A-B-C+E 

XDOT(3)=X(4) 

XDOT(4)=-P-Q+S 
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99  RETURN 

END 
C 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

SUBROUTINE  RK4 (T, X, XDOT, N, H, INDEX) 

IMPLICIT  REAL*8 (A-H,L-Z) 

INTEGER  INDEX, I 

DIMENSION  X(4) , XDOT (4) , SAVED (4) ,SAVEX(4) 
C 
C 

INDEX=INDEX+1 

GO  TO  (1,2,3,4), INDEX 

1  DO  10  1=1, N 
SAVEX (I) =X (I) 
SAVED ( I ) =XDOT ( I ) 

10     X(I)=SAVEX (I) +.5D0*H*XDOT (I) 

T=T+.5D0*H 

RETURN 
C 
C 

2  DO  20  1=1,  N 
SAVED(I)=SAVED(I)+2.D0*XDOT (I) 

20     X(I)=SAVEX(I) +.5D0*H*XDOT(I) 

RETURN 
C 
C 

3  DO  30  1  =  1, N 
SAVED(I)=SAVED (I) +2 .D0*XDOT ( I ) 

30     X(I)=SAVEX(I)+H*XDOT(I) 

T=T+.5D0*H 

RETURN 
C 
C 

4  DO  40  1=1, N 

4  0     X(I)=SAVEX(I)+H/6.D0* (SAVED (I) +XDOT(I) ) 

INDEX=0 

RETURN 

END 
C 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

SUBROUTINE  ORBDAT (ENERGY, ANGM, PI , a, e, APOGEE, PERIGE, PERIOD) 

IMPLICIT  REAL*8 (A-H, L-Z) 
C 
C 

e= (1+ (2*ENERGY*ANGM*ANGM/ (3 . 98601208133D+14* 
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*3. 98 6012  08 133D+14) ) ) **0.5 
a=(-3.9  8  6D+14/ (2*ENERGY) ) /1000 
APOGEE=a* (1+e) 
PERIGE=a* (1-e) 
PERIOD=( (2*PI) / (3.98  6D+05) **0.5)*(a**1.5) 


RETURN 

END 
C 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

SUBROUTINE  ONOFF (R0,X,TMAX, SPI, GO, TH, E0 , GAMMAR, ALPHAR, D) 

IMPLICIT  REAL*8 (A-H, L-Z) 

DIMENSION  X(4),XDOT(4) 
C 
C 

DELTAR=R0-1000 

V=(  (X(2)*X(2)  )  +  (X(l)*X(4)  )**2)**0.5 

MU=3.986012D+14 

M=2000  0 

E=M* ( (V*V) 12- (MU/X(1)  )  ) 
C 

C  PRINT*, 'MU=' ,MU 
C  PRINT*, 'E0  =  ' ,E0 
C 

IF  (TH  .EQ.  TMAX)  GO  TO  99 
C 

IF  (X(l)   .LE.  DELTAR)  THEN 
C 

IF  (E  .LT.  E0)  THEN 
TH=TMAX 


ELSE 
END  IF 


ELSE 
END  IF 
TH=0D+0 


TH=0D+0 


C 
C 

99  IF  (E  .GE.  E0)  THEN 

TH=0D+0 
ENDIF 

C 
C 

100  RETURN 
END 
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c 

c 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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